version 12.1
set more off
program drop _all
  quietly log
  local logon = r(status)
  if "`logon'" == "on" { 
	log close 
	}
log using bds2020psrm-postest.log, text replace


/*		********************************************************************	*/
/*     	File Name:		bds2020psrm-postest.do									*/
/*     	Date:   		March 03, 2020											*/
/*      Author: 		Frederick J. Boehmke									*/
/*      Purpose:		Do post-estimation analysis of the results from 		*/
/*						bds2020psrm.do. This creates graphs interpreting the 	*/
/*						war of attrition estimates for Figures 3-6 in Boehmke, 	*/
/*						Dion, and Shipan (2020, PSRM). 							*/
/*      Input File:		bds2020psrm-filibuster.dta								*/
/*						bds2020psrm-crises.xls									*/
/*						bds2020psrm-TXMY.ster									*/
/*      Output File:	bds2020psrm-postest.log									*/
/*						bds2020psrm-figure03.gph								*/
/*						bds2020psrm-figure04.gph								*/
/*						bds2020psrm-figure05.gph								*/
/*						bds2020psrm-figure06.gph								*/
/*		********************************************************************	*/

	/**************************************/
	/**************************************/
	/* Interpret the filibuster results.  */
	/**************************************/
	/**************************************/
	
		/* Figure 3: War of attrition model with no covariates. */

clear
set seed 540

		/* Open the saved estimates. */

estimates use bds2020psrm-T1M3

		/* Create a data set with t from 0 to 40. */

set obs 320

generat t = _n/8
  
		/* Create 10,000 copies. */

  expand 10000
  
  sort t
  
	bysort t: generat sim = _n
  
		/* Draw the parameters of the model from estimates. */

  matrix B = e(b)
  matrix C = e(V)
  
  drawnorm W_cons SS_cons logit_cons ln_p_W ln_p_SS, mean(B) cov(C)
  
  generat lambda_W = exp(W_cons)
  generat lambda_SS = exp(SS_cons)
  generat p_W = exp(ln_p_W)
  generat p_SS = exp(ln_p_SS)
  generat tstar = e(tstarmax)
  
		/* Calculate the estimated hazards for strong-strong */
		/* contests and contests involving at least one weak player. */
  
  generat haz_W = (lambda_W*p_W*(lambda_W*t)^(p_W-1))*((exp(-(lambda_W*t)^p_W))/(exp(-(lambda_W*t)^p_W) ///
	- exp(-(lambda_W*tstar)^p_W))) if t < tstar
  generat haz_SS = (lambda_SS*p_SS*(lambda_SS*t)^(p_SS-1))

	sort t

		/* Get percentiles of the draws for each value of t to construct CIs. */
	
  by t: egen haz_W_hi = pctile(haz_W), p(97.5)
  by t: egen haz_W_lo = pctile(haz_W), p(2.5)

  by t: egen haz_SS_hi = pctile(haz_SS), p(97.5)
  by t: egen haz_SS_lo = pctile(haz_SS), p(2.5)

		/* Collapse the data to have one observation per value of t. */
  
  collapse (mean) haz_W* haz_SS* tstar, by(t)

		/* Capture values we want to use in the graph. */
  
local lambda_W = exp([Weak]_b[_cons])
local p_W = exp([ln_p_W]_b[_cons])

local lambda_SS = exp([StrongStrong]_b[_cons])
local p_SS = exp([ln_p_SS]_b[_cons])

local prSS : display %3.2f (1+exp(-[logit_pi]_b[_cons]))^(-1)

		/* Now graph the results. */

twoway rspike haz_SS_lo haz_SS_hi t, scheme(s1mono)  /// 
		lcolor(gs10) lwidth(vthin) /// 
    || 	function y=`lambda_SS'*`p_SS'*((`lambda_SS'*x)^(`p_SS'-1)), range(1 40)  /// 
		lcolor(gs0) lwidth(medthick) lpattern(solid) /// 
    || 	rspike haz_W_lo haz_W_hi t if t < tstar-1, /// 
		lcolor(gs10) lwidth(vthin)  /// 
    || 	function y=(`lambda_W'*`p_W'*(`lambda_W'*x)^(`p_W'-1))*((exp(-(`lambda_W'*x)^`p_W'))/(exp(-(`lambda_W'*x)^`p_W')  /// 
			- exp(-(`lambda_W'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-1') /// 
		lcolor(gs0) lwidth(medthick) /// 
		ylabel(0(0.25)1.25, grid) ytitle(Hazard Rate) /// 
		xtitle(Length of Filibuster in Days) /// 
		xline(`e(tstarmax)', lcolor(gs10) lwidth(medthin)) /// 
		text(1.10 9  "Hazard for contests with at least one weak player", placement(e)) /// 
		text(0.22 20 "Hazard for Strong-Strong contests", placement(s)) /// 
		text(0.60 40 "Pr(Both Players Strong)= `prSS'", placement(nw)) /// 
		graphregion(margin(small)) /// 
		legend(off) /// 
		saving(bds2020psrm-figure03, replace)
	
		/* Figure 4: Now with covariates. */
		/* This plot does not have CIs because it was too hard to represent. */
	
estimates use bds2020psrm-T1M4

		/* Save the parameter values into locals. */

  local p_W = exp([ln_p_W]_b[_cons])
  
  local lambda_W1 = exp([Weak]_b[_cons])
  local lambda_W2 = exp([Weak]_b[_cons]+[Weak]_b[unclearbinder])  
  local lambda_W3 = exp([Weak]_b[_cons]+[Weak]_b[yesbinder])
  local lambda_W4 = exp([Weak]_b[_cons]+[Weak]_b[post1975])
  
  local lambda_SS = exp([StrongStrong]_b[_cons])
  local p_SS = exp([ln_p_SS]_b[_cons])

  local prSS1 : display %3.2f 1 - invlogit([logit_pi]_b[_cons])
  local prSS2 : display %3.2f 1 - invlogit([logit_pi]_b[_cons]+[logit_pi]_b[unclearbinder])  
  local prSS3 : display %3.2f 1 - invlogit([logit_pi]_b[_cons]+[logit_pi]_b[yesbinder])
  local prSS4 : display %3.2f 1 - invlogit([logit_pi]_b[_cons]+[logit_pi]_b[post1975])
  
		/* Now graph. */
  
  twoway function y=`lambda_SS'*`p_SS'*((`lambda_SS'*x)^(`p_SS'-1)), range(1 40) scheme(s1mono)  /// 
		lcolor(gs0) lwidth(thick) lpattern(solid) /// 
    || 	function y=(`lambda_W1'*`p_W'*(`lambda_W1'*x)^(`p_W'-1))*((exp(-(`lambda_W1'*x)^`p_W'))/(exp(-(`lambda_W1'*x)^`p_W') ///
			- exp(-(`lambda_W1'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-1') /// 
		lcolor(gs0) lwidth(thick) lpattern(solid) /// 
    || 	function y=(`lambda_W2'*`p_W'*(`lambda_W2'*x)^(`p_W'-1))*((exp(-(`lambda_W2'*x)^`p_W'))/(exp(-(`lambda_W2'*x)^`p_W') ///
			- exp(-(`lambda_W2'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-1') /// 
		lcolor(gs6) lwidth(thick) lpattern(solid) /// 
    || 	function y=(`lambda_W3'*`p_W'*(`lambda_W3'*x)^(`p_W'-1))*((exp(-(`lambda_W3'*x)^`p_W'))/(exp(-(`lambda_W3'*x)^`p_W') ///
			- exp(-(`lambda_W3'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-1') /// 
		lcolor(gs10) lwidth(thick) lpattern(solid) /// 
    || 	function y=(`lambda_W4'*`p_W'*(`lambda_W4'*x)^(`p_W'-1))*((exp(-(`lambda_W4'*x)^`p_W'))/(exp(-(`lambda_W4'*x)^`p_W') ///
			- exp(-(`lambda_W4'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-1') /// 
		lcolor(gs0) lwidth(thick) lpattern(longdash) /// 
		ylabel(#5, grid) ytitle(Hazard Rate) /// 
		xlabel(0 10 `e(tstarmax)' 20 30 40) /// 
		xtitle(Length of Filibuster in Days) /// 
		xline(`e(tstarmax)', lcolor(gs10) lwidth(medthin)) /// 
		text(2 -0.5 "Hazards for contests with" "at least one weak player", placement(se)) /// 
		text(0.1 40 "Hazard for Strong-Strong contests", placement(nw)) /// 
		text(1.0  14.1 "Important Issues, Pr(W|Imp) = `prSS3'", placement(se) color(gs10)) /// 
		text(1.4  14.1 "Unclear Importance Issues, Pr(W|Uncl) = `prSS2'", placement(e) color(gs6)) /// 
		text(1.7  14.1 "Unimportant Issues, Pr(W|Unimp) = `prSS1'", placement(e) color(gs0)) /// 
		text(2  14.1 "1975 Cloture Reform, Pr(W|CR) = `prSS4'", placement(se) color(gs0)) /// 
		graphregion(margin(tiny)) /// 
		legend(off) /// 
		saving(bds2020psrm-figure04, replace)

	
	/********************************************************************************/
	/* Do the PP plot for goodness of fit for filibusters (not reported in paper). 	*/
	/********************************************************************************/

		/* Model with no covariates. */
	
use bds2020psrm-filibuster, clear

		/* Calculate the predicted cdf for our estimator. */
  
estimates use bds2020psrm-T1M3

		/* Capture the parameters from the results. */

  predict xb_W, xb equation(Weak)
  predict xb_SS, xb equation(StrongStrong)
  predict xb_prSS, xb equation(logit_pi)

  generat lambda_W = exp(xb_W)
  generat lambda_SS = exp(xb_SS)
  
  local p_W = exp([ln_p_W]_b[_cons])
  local p_SS = exp([ln_p_SS]_b[_cons])
  local tstarmax = e(tstarmax)

		/* Calculate the cdfs for the two types of contests and then */
		/* average across them according to the mixture probability. */
  
  generat F_W = (1 - exp(-(lambda_W*days)^`p_W'))/(1 - exp(-(lambda_W*`tstarmax')^`p_W'))
  generat F_SS = 1 - exp(-(lambda_SS*days)^`p_SS')
  generat prSS = invlogit(xb_prSS)

  replace F_W = 1 if days > `tstarmax'
  
  generat F_comb = (1 - prSS)*F_W + prSS*F_SS
  
		/* Now do the Weibull. */
  
estimates use bds2020psrm-T1M1

  predict xb_wbl, xb
  generat lambda_wbl = exp(xb_wbl)

  local p_wbl = exp([ln_p]_b[_cons])
  
  generat F_wbl = 1 - exp(-(lambda_wbl*days)^`p_wbl')
    
		/* Calculate the empirical distribution of the sample (KM). */
  
  egen num_days = count(days), by(days)
  
  egen f_nonp = pc(num_days), prop
  
  collapse (sum) f_nonp (mean) F_*, by(days)
  
  sort days
  
  generat F_nonp = sum(f_nonp)
  
	correlate F_nonp F_comb F_wbl
 
		/* Graph the results. */

twoway function y=x, scheme(s1mono) /// 
		lcolor(black) lwidth(thin) lpattern(solid) /// 
  || 	connected F_comb F_wbl F_nonp, sort /// 
		lcolor(black gs6) lwidth(medthick medthick) lpattern(solid longdash)  /// 
		mlcolor(black gs6) mfcolor(black white) msize(small small) msymbol(D D) /// 
		xlabel(0(0.2)1, grid) ylabel(0(0.2)1, grid) /// 
		ytitle("Predicted") xtitle(Observed) /// 
		text(0.60 .60 "War of" "Attrition", place(se)) /// 
		text(0.47 0.82 "Weibull", place(se) color(gs6)) /// 
		legend(off) /// 
		xsize(5) ysize(5) 

	
	/**************************************/
	/**************************************/
	/* Interpret the crises results.      */
	/**************************************/
	/**************************************/

		/* Figure 5: The estimated hazard for crises. */
	
clear
set seed 540

		/* Open the saved estimates. */

estimates use bds2020psrm-T2M3

set obs 150

  generat t = _n/(200/300)
  
		/* Create 10,000 copies. */

  expand 10000
  
  sort t
  
	bysort t: generat sim = _n
  
  matrix B = e(b)
  matrix C = e(V)
  
		/* Draw the parameters of the model from estimates. */

  drawnorm W_cons SS_cons logit_cons ln_p_W ln_p_SS, mean(B) cov(C)
  
  generat lambda_W = exp(W_cons)
  generat lambda_SS = exp(SS_cons)
  generat p_W = exp(ln_p_W)
  generat p_SS = exp(ln_p_SS)
  generat tstar = e(tstarmax)
  
		/* Calculate the estimated hazards for strong-strong */
		/* contests and contests involving at least one weak player. */
  
  generat haz_W = (lambda_W*p_W*(lambda_W*t)^(p_W-1))*((exp(-(lambda_W*t)^p_W))/(exp(-(lambda_W*t)^p_W) ///
	- exp(-(lambda_W*tstar)^p_W))) if t < tstar
  generat haz_SS = (lambda_SS*p_SS*(lambda_SS*t)^(p_SS-1))

		/* Get percentiles of the draws for each value of t to construct CIs. */
	
	sort t

  by t: egen haz_W_hi = pctile(haz_W), p(97.5)
  by t: egen haz_W_lo = pctile(haz_W), p(2.5)

  by t: egen haz_SS_hi = pctile(haz_SS), p(97.5)
  by t: egen haz_SS_lo = pctile(haz_SS), p(2.5)

		/* Collapse the data to have one observation per value of t. */
  
  collapse (mean) haz_W* haz_SS* tstar, by(t)
  
		/* Capture values we want to use in the graph. */
  
generat haz_W_hi_cut = min( haz_W_hi, 0.4)

local prSS : display %3.2f (1+exp(-[logit_pi]_b[_cons]))^(-1)
local ll : display %3.2f e(ll)

local lambda_W = exp([Weak]_b[_cons])
local p_W = exp([ln_p_W]_b[_cons])

local lambda_SS = exp([StrongStrong]_b[_cons])
local p_SS = exp([ln_p_SS]_b[_cons])

		/* Now graph the results. */

twoway rspike haz_SS_lo haz_SS_hi t, scheme(s1mono)  /// 
		lcolor(gs10) lwidth(vthin) /// 
  || 	function y=`lambda_SS'*`p_SS'*((`lambda_SS'*x)^(`p_SS'-1)), range(1 300) /// 
		lcolor(gs0) lwidth(medthick) lpattern(solid) /// 
  || 	rspike haz_W_lo haz_W_hi_cut t if t <= tstar-5, /// 
		lcolor(gs10) lwidth(vthin)  /// 
  || 	function y=(`lambda_W'*`p_W'*(`lambda_W'*x)^(`p_W'-1))*((exp(-(`lambda_W'*x)^`p_W'))/(exp(-(`lambda_W'*x)^`p_W') - exp(-(`lambda_W'*e(tstarmax))^`p_W'))),  /// 
		range(1 `=`e(tstarmax)'-3') /// 
		lcolor(gs0) lwidth(medthick) /// 
		ylabel(#5, grid) ytitle(Hazard Rate) /// 
		xtitle(Length of Crisis in Days) /// 
		xline(`e(tstarmax)', lcolor(gs10) lwidth(medthin)) /// 
		text(0.35 44 "Hazard for contests with" "at least one weak player", placement(e)) /// 
		text(0.01 150 "Hazard for Strong-Strong contests", placement(n)) /// 
		text(0.175 300 "Pr(Both Players Strong)= `prSS'", placement(w)) /// 
		graphregion(margin(small)) /// 
		legend(off) /// 
		saving(bds2020psrm-figure05, replace)
		
	
	/************************************************************/
	/* Figure 6: The PP plot for goodness of fit for crises. 	*/
	/************************************************************/


		/* Model with no covariates. */
	
import excel using bds2020psrm-crises.xls, clear firstrow case(lower)

	generat logrelcap = ln(relcap)
	
		/* Calculate the predicted cdf for our estimator. */
  
estimates use bds2020psrm-T2M3

		/* Capture the parameters from the results. */

  predict xb_W, xb equation(Weak)
  predict xb_SS, xb equation(StrongStrong)
  predict xb_prSS, xb equation(logit_pi)

  generat lambda_W = exp(xb_W)
  generat lambda_SS = exp(xb_SS)
  
  local p_W = exp([ln_p_W]_b[_cons])
  local p_SS = exp([ln_p_SS]_b[_cons])
  local tstarmax = e(tstarmax)

		/* Calculate the cdfs for the two types of contests and then */
		/* average across them according to the mixture probability. */
  
  generat F_W = (1 - exp(-(lambda_W*trgterra)^`p_W'))/(1 - exp(-(lambda_W*`tstarmax')^`p_W'))
  generat F_SS = 1 - exp(-(lambda_SS*trgterra)^`p_SS')
  generat prSS = invlogit(xb_prSS)

  replace F_W = 1 if trgterra > `tstarmax'
  
  generat F_comb = (1 - prSS)*F_W + prSS*F_SS
  
		/* Now do the Weibull. */
  
estimates use bds2020psrm-T2M1

  predict xb_wbl, xb
  generat lambda_wbl = exp(xb_wbl)

  local p_wbl = exp([ln_p]_b[_cons])
  
  generat F_wbl = 1 - exp(-(lambda_wbl*trgterra)^`p_wbl')
    
		/* Now do the empirical distribution of the sample (KM). */
  
  egen num_trgterra = count(trgterra), by(trgterra)
  
  egen f_nonp = pc(num_trgterra), prop
  
  collapse (sum) f_nonp (mean) F_*, by(trgterra)
  
  sort trgterra
  
  generat F_nonp = sum(f_nonp)
  
	correlate F_nonp F_comb F_wbl
 
		/* Graph the results. */

twoway function y=x, scheme(s1mono) /// 
		lcolor(black) lwidth(thin) lpattern(solid) /// 
  || 	connected F_comb F_wbl F_nonp, sort /// 
		lcolor(black gs6) lwidth(medthick medthick) lpattern(solid longdash)  /// 
		mlcolor(black gs6) mfcolor(black white) msize(vsmall vsmall) msymbol(D D) /// 
		xlabel(0(0.2)1, grid) ylabel(0(0.2)1, grid) /// 
		ytitle("Predicted") xtitle(Observed) /// 
		text(0.56 .58 "War of" "Attrition", place(se)) /// 
		text(0.28 0.64 "Weibull", place(e) color(gs6)) /// 
		legend(off) /// 
		xsize(5) ysize(5) /// 
		saving(bds2020psrm-figure06, replace)
	
log close
clear
exit, STATA
